
#include <util/common_macros.h>
#include <util/math_macros.h>

#include "../../imp_types.h"
#include "../../matrix/matrix.h"
/**
 * \brief Calculates the integral image.
 * 
 * Calculates the integral (Summed Table Area) of an image.<br>
 * <br>
 * This function calculates integral image as following:<br>
 * 
 * \f[
 * IntImg(X,Y) = \sum_{x < X, y < Y} Img(x,y)
 * \f]
 * 
 * @param Img input matrix
 * @param IntImg the matrix that will be the integral image. If \c *IntImg is equal to NULL, \c *IntImg will be allocated in this function.
 * 
 */
void imp_image_integral(ImpMatrix * Img, ImpMatrix ** IntImg) {

	real * img;
	real * int_img;

	int r, c, _r, i;
	int rows, columns;
	real sum;

	if (*IntImg == NULL) (*IntImg) = imp_matrix_new(Img->rows, Img->columns);

	img = Img->data;
	rows = Img->rows;
	columns = Img->columns;

	int_img = (*IntImg)->data;

	sum = 0;

	for (i = 0; i < columns; i++) {

		sum = sum + img[i];
		int_img[i] = sum;

	}
	for (r = 1, _r = r * columns; r < rows; r++, _r += columns) {

		sum = 0;

		for (c = 0; c < columns; c++) {
			i = _r + c;

			sum = sum + img[i];
			int_img[i] = int_img[i - columns] + sum;

		}
	}
}
/**
 * \brief Calculates the integral image.
 * 
 * Calculates the integral (Summed Table Area) of an image.<br>
 * <br>
 * This function calculates integral image as following:<br>
 * 
 * \f[
 * IntImg(X,Y) = \sum_{x < X, y < Y} Img(x,y)
 * \f]
 * 
 * \f[
 * SqIntImg(X,Y) = \sum_{x < X, y < Y} Img(x,y)^{2}
 * \f]
 * 
 * @param Img input matrix
 * @param IntImg the matrix that will be the integral image. If \c *IntImg is equal to NULL, \c *IntImg will be allocated in this function.
 * @param SqIntImg integral image of squared pixels. If \c *SqIntImg is equal to NULL, \c *SqIntImg will be allocated in this function.
 * 
 */
void imp_image_integral2(ImpMatrix * Img, ImpMatrix ** IntImg, ImpMatrix ** SqIntImg) {

	real * img;
	real * int_img;
	real * sq_int_img;

	int r, c, _r, i;
	int rows, columns;
	real sum, sqsum;

	if (*IntImg == NULL) (*IntImg) = imp_matrix_new(Img->rows, Img->columns);
	if (*SqIntImg == NULL) (*SqIntImg) = imp_matrix_new(Img->rows, Img->columns);

	img = Img->data;
	rows = Img->rows;
	columns = Img->columns;

	int_img = (*IntImg)->data;
	sq_int_img = (*SqIntImg)->data;

	sum = 0;
	sqsum = 0;

	for (i = 0; i < columns; i++) {

		sum = sum + img[i];
		int_img[i] = sum;

		sqsum = sqsum + (img[i] * img[i]);
		sq_int_img[i] = sqsum;

	}
	for (r = 1, _r = r * columns; r < rows; r++, _r += columns) {
		sum = 0;
		sqsum = 0;
		for (c = 0; c < columns; c++) {
			i = _r + c;

			sum = sum + img[i];
			int_img[i] = int_img[i - columns] + sum;

			sqsum = sqsum + (img[i] * img[i]);
			sq_int_img[i] = sq_int_img[i - columns] + sqsum;
		}
	}
}
/**
 * \brief Calculates the area of a rectangle of an integral image.
 * 
 * Calculates the area of a rectangle of an integral image.
 * 
 * @param Integral Integral of the image.
 * @param p1 top-left point coordinate.
 * @param p2 downl-right point coordinate.
 * @return area of rectangle.
 */
real imp_image_integral_get_area(ImpMatrix * Integral, Ref2 p1, Ref2 p2) {
	real area;

	p1.x = CLAMP(p1.x, 0, Integral->columns - 1);
	p1.y = CLAMP(p1.y, 0, Integral->rows - 1);

	area = Integral->data[ imp_matrix_idx(Integral, p2.x, p2.y) ];
	area = area - Integral->data[ imp_matrix_idx(Integral, p1.x, p2.y) ];
	area = area - Integral->data[ imp_matrix_idx(Integral, p2.x, p1.y) ];
	area = area + Integral->data[ imp_matrix_idx(Integral, p1.x, p1.y) ];
	return area;
}
/**
 * \brief Calculates the area of a rectangle of an integral image.
 * 
 * Calculates the area of a rectangle of an integral image.
 * 
 * @param Integral Integral of the image.
 * @param rect a rectangle in the integral image.
 * @return area of rectangle.
 */
real imp_image_integral_get_area_from_rect(ImpMatrix * Integral, Rect rect) {
	return imp_image_integral_get_area(Integral, rect.p1, rect.p2);
}
